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Abstract 

The computation of the ground state (i.e. the eigenvector related to the small- 
est eigenvalue) is an important task in the simulation of quantum many-body 
systems. As the dimension of the underlying vector space grows exponentially 
in the number of particles, one has to consider appropriate subsets promis- 
ing both convenient approximation properties and efficient computations. The 
variational ansatz for this numerical approach leads to the minimization of the 
Rayleigh quotient. The Alternating Least Squares technique is then applied to 
break down the eigenvector computation to problems of appropriate size, which 
can be solved by classical methods. Efficient computations require fast compu- 
tation of the matrix- vector product and of the inner product of two decomposed 
vectors. To this end, both appropriate representations of vectors and efficient 
contraction schemes are needed. 

Here approaches from many-body quantum physics for one-dimensional and 
two-dimensional systems (Matrix Product States and Projected Entangled Pair 
States) are treated mathematically in terms of tensors. We give the definition of 
these concepts, bring some results concerning uniqueness and numerical stability 
and show how computations can be executed efficiently within these concepts. 
Based on this overview we present some modifications and generalizations of 
these concepts and show that they still allow efficient computations such as ap- 
plicable contraction schemes. In this context we consider the minimization of 
the Rayleigh quotient in terms of the PARAFAC (CP) formalism, where we also 
allow different tensor partitions. This approach makes use of efficient contrac- 
tion schemes for the calculation of inner products in a way that can easily be 
extended to the MPS format but also to higher dimensional problems. 

Keywords: Quantum many-body systems, Density Matrix Renormalization 
Group (dmrg), Matrix Product States (mps), Tensor Trains, Projected 
Entangled-Pair States, Canonical Decomposition (candecomp or parafac) 
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1. Introduction 

Computations with tensors are getting increasingly important in high di- 
mensional problems. In particular in quantum many-body physics, a typical 
problem amounts to finding an accurate approximation to the smallest eigen- 
value (i.e., the ground-state energy) of a hermitian matrix (representing the 
Hamiltonian) that is larger than one can store even on a powerful computer. To 
this end, in quantum physics techniques like Matrix Product States (mps) or 
Projected Entangled Pair States (peps) have been developed for representing 
vectors, viz. eigenstates of quantum systems efficiently. In mathematics, besides 
the Tucker decomposition and the canonical decomposition, concepts like Ten- 
sor Trains (tt) were introduced. The examples of MPS or PEPS (in physics) and 
tt (in mathematics) express a common interest in powerful numerical meth- 
ods specifically designed for coping with high-dimensional tensor networks. - 
Unifying variational approaches to ground-state calculations [l| in a common 
framework of tensor approximations will be highly useful, in particular in view 
of optimizing numerical algorithms 0, S 0, Q • Here it is the goal to cast some 
of the recent developments in mathematical physics into such a common frame 
expressed in the terminology of multilinear algebra. Moreover we present nu- 
merical results on the Ising-type Hamiltonian underscoring the wealth and the 
potential of such an approach. 

In this paper, we address one-dimensional and two-dimensional methods in a 
unified frame related to MPS and peps. We will introduce some generalization of 
MPS and the canonical decomposition. Furthermore, we give a short description 
of tensor-decomposition methods for 2D problems. 

Scope and Organization 

The paper is organized as follows: Sections [2] [3] and [4] contain an overview of 
already-known concepts and describe them in the multilinear algebra language: 
in Section [2] we introduce the physical background of the problem setting and 
define the matrices involved, in Section |3] we present representation schemes for 
states in physically motivated ID and 2D arrays and we show how computations 
can be performed efficiently and Section@]finally fixes some basics and notations 
for tensors and tensor-decomposition schemes. 

In Section [5] we present new ideas of how to generalize these basic concepts, 
how to execute calculations efficiently and how to apply them to the ground- 
state approximation problem. First numerical results will show the benefit of 
these newly developed concepts. 

2. Physical Model Systems 

Consider vectors a; in a complex Hilbcrt space T-L representing states of (pure) 
quantum systems. The differential equation x = —iHx (Schrodinger's equation 
of motion) then governs quantum dynamics (neglecting relaxation) with the 
Hamiltonian H being the generator of unitary time evolution. The Hamiltonian 
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captures the energies of the constituent subsystems (e.g. spins) as well as the 
interaction energies between coupled subsystems. 

For instance, a linear chain of five spins coupled by nearest neighbor interac- 



tions can be depicted as in Figure 1(a) In the case of open boundary conditions 
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(a) ID system of 5 spins with nearest-neighbor interaction and periodic boundary condi- 
tions. 
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(b) Representation of the Hamiltonian related to 
the physical system illustrated by Figure [l (a) | 

Figure 1: Example of a linear physical system (a) and the related Hamiltonian (b). 

(OBC) there is no coupling interaction between particle 1 and 5, while for pe- 
riodic boundary conditions (PBC) there is a non vanishing coupling interaction 
between 1 and 5. 



2.1. Hamiltonian Representations 

For spin | particles such as electrons or protons, the spin angular momentum 
operator describing their internal degree of freedom (i.e. spin- up and spin-down) 
is usually expressed in terms of the Pauli matrices 



P, 



1\ p 

1 or y 



o 







and P z 



1 

-1 



(1) 



Being traceless and Hermitian, {P x , P y , P z } forms a basis of the Lie algebra 
su(2), while by appending the identity matrix I one obtains a basis of the Lie 
algebra u(2). 

Now, spin Hamiltonians are built by summing M terms, each of them rep- 
resenting a physical (inter) action. These terms are themselves tensor products 
of Pauli matrices or identities 

M M 

H = J2 MQ[ k) ® Q ( 2 k) ® • • • ® Qf ] ) = E H (k) > ( 2 ) 



k=l 



= :H<. k > 



k=l 
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(k) 

where Qj ' can be P x , P y , P z or I. 

In each summand most of the arc /: local terms have just one non- 
trivial tensor factor, while pair interactions have two of them. Higher m-body 
interactions (with m > 2) usually do not occur as physical primitives, but could 
be represented likewise by m Pauli matrices in the tensor product representing 
the ?7i-order interaction terrr0. 

For instance, in the Ising (ZZ) model [3] for the ID chain with p spins and 
open boundary conditions, the spin Hamiltonian takes the form 

p-i 

fcl p (3) 

+ A /® (fc_1) ® (P,)* (8 /®(P- fc ) , 
fe=i 

where the index fc denotes the position in the spin chain and the real num- 
ber A describes the ratio of the strengths of the magnetic field and the pair 
interactions. For simplicity, we will henceforth drop the tensor powers of the 
identity and tacitly assume appropriate embedding. Then a Hamiltonian for an 
open-boundary ID Heisenberg (XY) model [1, Q reads 

p-i 

H = J x ■ J (8 (P x ) k <E) (P x ) k +1 ®I + JyI® (Py)k <8 (Py)k+1 8> I 



k=l 



(4) 



®{p x ) k ®I 



k=l 



with real constants J x , J y and A. Models with all coefficients being different are 
called anisotropic. 

Being a sum ((2|) of Kronecker products of structured 2x2 matrices many 
Hamiltonians have special properties: they can be multilevel-circulant ( 1^, lH) 



or skew-circulant, diagonal or persymmctric (|12|), which can be exploited to 
derive properties of the respective eigenvalues and eigenvectors. 



2.2. Computation of Ground States: Physical Background 

A key to understand some of the motivating guidelines lies in the somewhat 
striking fact that quantum dynamical systems typically evolve in a way that 
looks non-generic from a mathematical point of view. Yet the very structure 
of quantum dynamics paves the way to tailored parameterizations based on 
tensor compression that arc efficient in the sense of scaling only polynomially 
in physical system size. Some of these motivating guidelines established in 
quantum physics may be sketched as follows: Composing a quantum system 



1 For further details, a reader wishing to approach quantum physics from linear and multi- 
linear algebra may refer to [f|. 
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from its components takes a joint Hilbert space that is the tensor product of 
the individual Hilbert spaces. Likewise a linear operator on the joint Hilbert 
space can be composed by taking sums (or weighted linear combinations) of 
tensor factors (like in Eqn.[2]). Clearly, in general a linear combination of tensor 
products does not take the form of a tensor product itself. Thus a quantum 
state space grows exponentially with the number of constituents in contrast to 
a classical configuration space just growing linearly. 

However, correlating quantum interactions typically become smaller and 
smaller with increasing distance between subsystems ('particles'): for instance, 
in Eqn. |3]only nearest-neighbor interactions had to be taken into account. On 
a general scale, this can be made precise in terms of area laws, where the cor- 
relations are quan tified by a measure termed 'entanglement entropy' of ground 
states [H, 14, lB[, see also the recent review in 16 1. Remarkably, this entropy 



of the reduced state of a subregion is not extensive: it typically grows with 
the boundary region ('area') between the subregion and its complement rather 
than with the volume of the subregion. In one-dimensional systems, a rigor- 
ous area law has recently been proven for a ll sy stems with a gap between the 
smallest and the second smallest eigenvalue [17[. Extending the results to two- 
dimensional lattice systems, however, currently requires stronger assumptions 
on the eigenvalue distribution [HI ■ 

Guided by these area laws, long-distance correlations may be neglected in 
the sense that eigenvectors (ground states) of physical systems are well approx- 
imated within truncated subsets, as has been quantitatively established, e.g., 
for MPS [lij]. Moreover MPS-approximations to ground states can provably be 
calculated efficiently [2(|. Along similar lines, consecutive partitionings have 
been exploited in unitary and tensor networks addressing ground states and dy- 
namics of large-scale quantum systems. Related techniques for truncating the 
Hilbert space to pertinent parameterized subsets not only include Matrix Prod- 
uct States (mps) [2l|, [22| of Density Matrix Renormalization Groups (dmrg) 



23l . l24f . but also pro jected entangled pair states (peps) [2a, [26j, weighted 
graph states (wgs) [27| , Multi-scale Entanglement Renormalization Approaches 
(mera) [Hj], string-bond states (sbs) [29| as well as combined methods fit [30l| . 

To conclude, the evolution of physical systems does not exploit the generic 
state space (with long-distance many-body interactions) , but proceeds via well- 
defined subspaces of short-distance and mainly pairwise interactions that can 
be parameterized by data-sparse formats which allow for tensor-contraction 
schemes. 

2.3. Computation of Ground States: Numerical Aspects 

The ground state energy of a physical system modeled by the Hamilto- 
nian H corresponds to the smallest eigenvalue of H, which is the minimum of 



the Rayleigh quotient [31 1 



min — ^ — . (5) 

xGH x H x 



As long as p, the number of particles, is not too large, any standard numerical 
method for computing eigenvalues and eigenvectors can be used. 
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But with increasing p, the size of the Hamiltonians grows like 2 P . Thus for 
p > 50 neither matrices nor even vectors of this size can be stored. So, similar 
to the description of the Hamiltonian @ we need a sparse approximate repre- 
sentation of eigenvectors. Assume that we have already chosen an appropriate 
subset UcH, the goal is to find approximations for the eigenvector in this set. 
Hence we consider the minimization of the Rayleigh quotient ([5]) only on the 
subset U: 

x H Hx , . 

mm — — . (6) 

xEU X H X 

An appropriate subset of vectors should allow for easy computation of Hx and 
of inner products y^x. Therefore we consider vector representations with a less 
number of coefficients where subsets of the indices can be grouped in partitions 
corresponding to the binary tensor structure of the Hamiltonian ([2]). Please 
note that, in general, the chosen subsets do not form linear subspaces. 



2-4- Alternating Least Squares 

An important tool for minimizing the Rayleigh quotient for an appropriate 
subset U is the Alternating Least Squares approach (als) , see 32J, |33| . Here, 
all subsets of the partitioning up to one are assumed to be fixed, and then the 
minimization is reduced to the remaining subset. As introductory example let 
us look at a set of vectors defined by 

X = X X ® X 2 <2> ■ • • ® X p = (X 1;il ■ X 2 :i 2 ■ ■ ■ Xp;i p )i u ...,ip = {Xi)i=0,- ,2»-l (7) 

with vectors Xi of length 2, and i = (ii, ... ,i p )2 the binary representation of 
i with ij G {0, 1}. Hence we envisage the vector s as a p-tensor. So in our 
example ([7]) we assume all subsets fixed up to x r , and then the minimization is 
simplified to 



x u Hx 



(sbi ® • • • <8> Xpf ( jr a k Q[ k) <g> • • • g) Q p k) ^j (xi <S> ■ ■ 
(xi ® ■ ■ ■ ® x p ) H (xi ® • • • <g> x p ) 

^ ] C^fc(«Cl Q\ ^l) ' ' ' {,^r Qr 3Cy ) ' ' ' {Xp Qp 



■ k=l 
mm 



E a k /3 k Qi k) ) x 



= mm ■ 



a standard eigenvalue problem in the effective Hamiltonian R r = J2k=i ^r^Qr^ ■ 
More generally, if we work with more complex representations such as Matrix 
Product States, the minimization of the Rayleigh quotient © will lead to 

. x»Hx . x»R r x r (8) 
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with an effective Hamiltonian related to the generalized eigenvalue problem 
in R r and N r . So far, Eq. [8] describes the general situation, the particular 
specification of the matrices R r and N r will be given when we consider different 
representation schemes for the eigenvector. Then, x r is set to the eigenvector 
with smallest eigenvalue. We can repeat this procedure step by step for all 
Xj to get approximations for the minimal eigenvalue of H. The main costs 
are caused by matrix-vector products Hx and inner products y^x plus the 
solution of the relatively small generalized eigenvalue problem (JHJ. It is therefore 
important to have efficient schemes for the evaluation of inner products of two 
vectors out of the chosen subset U. Wc emphasize that this approach (and 
adapted modifications) allows to overcome the curse of dimensionality as it is 
only polynomial in the maximum length of the small vectors Xj , in the number 
of such vectors (which can be upper bounded by p) and in the number of local 
terms M. 

The ansatz ([8]) may cause problems if the denominator matrix N r is singu- 
lar. In that case one would apply an orthogonal projection on the nonsingular 
subspacc of N r . 



3. Representations of States for Computing Ground States 

Obviously, in general the above simplistic approach based on a single tensor 
product cannot give good approximations to the eigenvector. Therefore, we 
have to find a clever combination of such terms. The choice naturally depends 
on the dimension and the neighborhood relation of the physical setting. So first 
we consider the ID linear setting, and in a following section we look at the 2D 
problem. 



3.1. ID Systems: Approximation by Matrix Product States 

The Matrix Product State (mps) formalism goes back to several sources: 
early ones are by Affleck, Kennedy, Lieb, and Tasaki [HI, Q including their 
revival by Fannes, Nachtergaele, and Werner 2l|,l22|, while a more recent treat- 
ment is due to Vidal [35|. The application to the eigenvalue problem was dis- 
cussed by Delgado et al. [36| . A seemingly independent line of thought resorts 
to the fact that Density Matrix Renormalization Group (DMRG) methods as 



developed by Wilson and White [37J, |38( have a natural interpretation as opti- 
mization in the class of mps states, see, e.g., Refs. [3, 0, Hsf. As has been 
mentioned alrea dy, g round states of gapped ID Hamiltonians arc faithfully rep- 
resented by MPS [19( , where the MPS-approximation can be computed efficiently 
"391 . 20j , the rationale being an area law [ItJ • 



3.1.1. Formalism and Computations 

For MPS small Dj x -matrices are used for describing vectors in a com- 
pact form. The advantage is due to the fact that D := max{Dj} (the bond 
dimension) has to grow only polynomially in the number of particles in order 
to approximate ground states with a given precision (flij). 
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In MPS, the vector components are given by 
x i = x iu ..., ip =tr(A^.At ) ---A^) 



(9) 



7 (^) n (i P ) 

u "l;mi,m2 L 

mi — 1 m„-l 



. n n v p - 

rri2 2;m2,m3 '^mp,??!! 



The matrix products lead to indices and summation over m2, ...m p , and the 
trace introduces mi. The upper (physical or given) indices identify which of 
the two possible matrices are used at each position, and thereby they determine 
the vector components. So, e.g., the last component is described by 

X2P-1 — •*•].,.. .,1 — '^2 i — / y tI l;mi,m 2 ' li 2;m2,m3 ' ' ' u p;m p ,m 1 ! 

mi ,. . . ,m p 

where we always choose the matrix index ij = 1. The additional summation 
indices are called ancilla indices. 

The above trace form is related to the periodic case. In the open boundary 
case, there is no connection between first and last particle, and therefore the 
index mi can be neglected. 

In this case we have D\ = D p+ i = 1 and thus the matrices at the ends are 
of size 1 x Z?2 and D p x 1 respectively. 

2-8!,...,%, — ^1 ■ p-1 P 



D 2 D p /-.qx 

Mi) .Ji2) (ip-i) {ip) 

l;l,m2 2;T7i2,m3 p— l;m p _i,m p p;m p ,l * 

mo— 1 m p — l 



E • • • ^2 a l-A,m 2 ' fl 2; 



By introducing the unit vectors = e^,...^ = ® • • • ® with unit 
vectors e^^. of length 2, another useful representation of the MPS vector is given 
by 



— E x *i ip e ti,.».i P — E tr (^i^ ■^■p , '') e ii,—,* 



(11) 



LI 

V V a (il) e 

ii,— ,i p mi,— ,m p 

mi,...,m p zi i p 

mi,- ,m p 

with length 2 vectors a r . rnr , : m r+1 where the two components are pairwise entries 
in the matrices and A^ at position (m r ,m r+ i): 



fl r;m r ,m r+1 : I (1) 



(0) 



, a 



r;m r ,m r +i 
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Uniqueness of MPS and Normal Forms 

In this section, we want to summarize some known results concerning the 
uniqueness of MPS. For further details, see, e.g., [4(|. Obviously the represen- 
tation of an MPS vector is not unique. So, for a vector with components 

*i = = tr (4*° • 4 2) 4-1° • (12) 

we can replace the matrices by 

Af ] -> M^A^Mj+t , A[ n) -> A ( ; i] M 2 , A^ -> M^A^ (13) 

with nonsingular matrices Mj 6 C Dj xDj ,j = 2, . . . , p. 

The absence of uniqueness also causes problems in the solution of the effective 
generalized eigenvalue problem, because the matrix N r in Eqn. (|5J) might be 
positive semidefinite, but singular. To avoid this problem, we switch from a 
given MPS representation to a representation based on unitary matrices. To this 
end, we combine the matrix pair A^ 1 ' for i r = 0, 1 to a rectangular matrix and 
compute the SVD: 

(o)-*-($)-<mh < 14 > 

where the Ur lr ^ are the left part of U r . Now we can replace at position r 
in the MPS vector the matrix pair 4 by the pair Ur ^ and multiply the 
remaining SVD factor A r V r from the left to the right neighbor pair 4+i 
without changing the vector: 

tr (4° • Af ] ■ ■ ■ AW ■ A^ ■■■■ 4-1° ■ A p p) ) ~ ► 
tr • 4 82) • • • C/W • (ArV^A^ ■ ■ ■ 4-1° ' 4 lp) ) ■ 

So we can start from the left, normalizing first A± , always moving the remain- 
fa ) 

ing SVD part to the right neighbor, until we reach A p " . During this procedure 
the MPS matrix pairs Aj, j = 1, . . .,p — 1 are replaced by parts of unitary 
matrices C/j lj , which fulfil the gauge condition 

Y,uf i)n ufi ] =I. (15) 

In the case of open boundary conditions, the right (unnormalized) matrices 
Ap p ^ are column vectors and thus Ap p ^ H Ap is only a scalar 7, which 
corresponds to the squared norm of the MPS vector: 

x *x = Yl (4 il) ---4 p) )-(4 il) ---4 ip) ) 
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ti,"- Ap 

E4 lp)H (•••E4 42)H (E4 ll)H 4 4l) ) 4 42) ---) 4^ 

i P \ J2 V ii / / 



E3 
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Thus, if a; has norm one, the gauge condition (|15|) is also fulhllcd for j = p. 

The same procedure can be applied in order to move the remaining SVD 
part to the left neighbor. To this end, we compute 

(4 0) Ai 1} ) = V r ■ (A r 0)-U r = (V r A r ) (u^ . (16) 

Similarly we can move from right to left and replace the matrix pairs Aj , 

j = p, . . . , 2 by the unitaries U^ 3 ^ until we reach A± . Now the gauge conditions 
take the form 

Y J Uf l) uf l)Yl =I . (17) 

i 3 

Analogously, for open boundary conditions the remaining left matrices A^ are 
row vectors and so J2i A[ l1 ^ A 1 ^ 1 ^ is simply a scalar, which is 1 for a norm 1 
vector. 

So far, for the normalization process only one matrix pair A^ was involved. 



Similar to the two-site DMRG approach [24[, it is also possible to consider 
the matrices related to two neighboring sites at once [41|. To this end, we 
consider the two matrix pairs A ( - 3) 6 C^* 13 * 1 and e C D t+ lXD '+ 2 . The 

four matrix products AjA^p^ <E £, D j xD i+2^ are now re-arranged in matrix 
notation and an SVD is carried out: 



(i) j \ a j+i A j+i) - 1 A w a(°) a [1) a (1) I ~ \ r/ (1) I 3 V 3+ 1 



If we sweep from left to right, we replace the matrices Aj by parts of unitary 

(i -) 

matrices ?/• , shift the remaining part to the right neighbor, i.e. 

and proceed with the adjacent sites j + 1 and j + 2. Accordingly, if we sweep 
from right to left, we replace the matrices A^j^ 1 ' by the unitaries vfe*" 1 ', shift 
the remaining part to site j, i.e. 
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and proceed with the index pair (j — l,j). 

There exist even stronger normalization conditions allowing representations 
which are unique up to permutations and degeneracies in the singular values, 
see, e.g. [HI, l4o| . The proof of the existence of such normal forms is based on 
the SVD of special matricizations of the vector to be represented, see 42. 43l|. 
In the SVD-TT algorithm (44[, the same technique is in use. In 43 1 we present 
normal forms for MPS which allow for expressing certain symmetry relations. 

The presented normalization techniques for MPS vectors have various ad- 
vantages. They introduce normal forms for MPS vectors which lead to better 
convergence properties. For the minimization of the Rayleigh quotient ([8]), the 
gauge conditions circumvent the problem of bad conditioned N r matrices in 
the denominator and therefore approve numerical stability. So far, the pre- 
sented normalization technique only changes the representation of the vector 
but does not change the overall vector. However, the SVD can also be used as 
a truncation technique. This could be interesting if we want to keep the matrix 
dimensions limited by some D = Z? max - As an example we mention the PEPS 
format [45| . where such an SVD-based reduction appears, compare Subsection 



Sum of MPS Vectors 

Unfortunately, the MPS formalism does not define a linear subspacc. Ac- 
cording to [4l[ the sum of two MPS vectors x and y, which are both in PBC 
form, can be formulated as 



<i *p 

E * 



E 



tr 



= E 

i 1 ,...,i p 



(ii) 



(u) 



■A 



(ip) 



(ip) 



(»l) 



»1j..-j»j> 



B[ n) ■ ■ ■ B ( p p) 



1 (ip) 



B 



eii, 



(ip) 



(18) 



eii, 



In the open boundary case, we can define the interior matrices C2, . . . , C p _i in 
the same way. For reasons of consistency, the C matrices at the boundary sites 
1 and p have to be specified to be vectors, i.e. 



a[ 1i \b[ 1i) 



(j(i P ) 



A. 



(ip) 



Hence, the sum of MPS vectors is again an MPS vector, but of larger size. In order 
to keep the sizes of the MPS matrices limited by a constant -D max , one could apply 
an SVD-based truncation of the resulting MPS formalism and consider only the 
-Dmax dominant terms. 
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Properties of MPS 

Every unit vector can be represented by an MPS with bond dimension D = 1: 
let j = (Ji, . . . ,j p ) be the binary form of j, then ej can be formulated as an 
MPS with lxl MPS matrices Ay r ' = 5 ir j r : 

«- t <^,^>. 

Zl,...,2p— 

In view of Eqn. (|18p this observation may be extended as follows: every sparse 
vector with at most D non zero entries can be written as an MPS with bond 
dimension D. 

Computation of Matrix- Vector Products 

To solve the minimization of the Rayleigh quotient ((SJ) in an ALS way as 
introduced in Subsection 12.41 we have to compute 



Vk = X = [oLkQ^ ® ■ ■ ■ ® Qf ] y J2 (°lim ll m s ® • • • ® 

mi , . . .,m p 

mi,. . .,m p 

= ^ afefbl.fejmi.ma ® •• • ® &p,fc;m p ,mi 

mi,...,m p 

and derive a sum of MPS vectors 

M M 

y = Hx = Y J H {k) x = Y j yk ■ (19) 



k=l k=l 



Therefore, we can compute all small vectors Qj k ' ^ a,j;m j ,m :j+1 in 0(2pD 2 M) 
operations. For reasons of efficiency there are also concepts to express the 
Hamiltonian in terms of a Matrix Product Operator (mpo) as defined in (|2"2"j) . 
This approach enables short representations for MPO x MPS products. 

Computation of Inner Products 

Furthermore, we have to compute inner products of two MPS vectors 

V V (b (il) b {ip) \-(a in) \ (20) 

/ 4 / , yi-MM U p;k p .k 1 J \ u l;m 1 ,m 2 u p;m p , mi J • ^ u ) 

To deal with these sums, it is essential to specify in which order the summations 
have to be conducted. In the next figures we display such an efficient ordering 

for the periodic case. To this end, we introduce the following notation. Each box 

(i ) 

in the following figures describes one factor, e.g., in these collections 

of sums. Little lines (legs) describe the connection of two such boxes via a 
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common index. Hence, the number of legs of a box is exactly the number of 



indices. So in this case, most of the boxes have three legs. Figure 2(a) shows 
the sum and the first terms, resp. boxes, with their connections. This also 
represents a Tensor Network. 

Now, in a first step we reorder the sums, and execute a contraction relative 
to index i\. This is given by the partial sum 



(»i) _ 

fe 2 a l;mi,m a — Ck 1 ,k 2 ;m lt m 2 ■ 



(21) 



as 



This eliminates two boxes, but leads to a new four leg tensor Ck 1 ,k 2 :m 1 ,m 2 
shown in |2(b)| and 3(a) Now we contract index fe, as shown in Figure [3(b)[ 



TT 



2*2, &2, ^3 








, ni2 , ni 3 
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ii,ti,fa 
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12, m 2 , TTI3 




i 3 , 7713 , m 4 
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(a) Computation of the inner product J 20 1 ) 
of two MPS vectors. 



(b) Contraction of the two MPS vectors con- 
cerning index i\. 



Figure 2: Contraction of the inner product of two MPS vectors. The dashed red line illustrates 
the index being contracted. 



leading to 3(c) In the following step we go from 3(c) to 3(d) by contracting i% 
and TO2, deriving at Figure |3(d) I Now we are in the same situation as in Figure 



3(a) and we can proceed exactly in the same way, until all sums are executed, 
resp. all indices have been contracted. 

The main costs depend on the size of the contracted index, e.g. 2 for con- 
tracting i r or D for contracting m r or k r , and on the size of the other indices 
that appear in the contracted boxes. Hence, e.g. the contraction in Figure 
|3(b)| costs D ■ 2D 3 = 2D i operations, and Figure |3(c)| costs 2D • D i = 2D 



operations. The total costs for the inner product is therefore less than 4£> 5 p, 
because contractions have to be done until all boxes are removed, that is 2p- 
times. Therefore, the costs for computing x^Hx based on (fT9|) arc less than 
AD 5 Mp. 

In the open boundary case, we start at the left or right side and therefore 



only contractions of boxes of length 3 can occur and thus Figure 3(c) only costs 
2D ■ D 2 = 2D 3 instead of 2D 5 in the periodic case. Thus, the overall costs for 
the inner product is less than 4D 3 p and for x H Hx less than AD 3 Mp. 



Minimization of the Rayleigh Quotient in Terms of MPS 

Now, we use the MPS formalism as vector ansatz to minimize the Rayleigh 
quotient (|6]). Therefore we start the ALS procedure, updating A p p by an im- 
proved estimate by solving the generalized effective eigenvalue problem. In 
addition, we replace A p p by a unitary pair via SVD, moving the remaining 
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r~"fe""i r 



i _l I L 

(a) After the ii-contraction, we get a four 
leg tensor 



i _l I L 

(b) Contraction concerning index £12 



*i,mi, i 2 ,m 2 ,fc. 



T 



T 



(c) Contraction concerning the indices £2 (d) After the contraction concerning , w.2 

and m2 and 12, we are in the same situation as in 

Figure [3(a)] 



Figure 3: Contraction of the inner product of two MPS vectors. 



SVD term to the left neighbor, and so on. A nice side effect here is that in the 
case of open boundary conditions the matrix N r is the identity because all the 
matrix pairs to index different from r are parts of unitary matrices and thus 
fulfil one of the gauge conditions ([T5|) or (fTTj) . Together with the fact that the 
MPS matrices at both ends are simply vectors we obtain 

x*x = ti(A^ ■ ■ ■ A^) ■ tr(A^ ■ ■ ■ A^) 

iir" Ap 

= Yl tr("(4 il) ---4 i " ) )®(4 ii) ---4 ip) )) 

i\ ,i p 

= tr((4 ll) ®4 n) )---(4 ip) ®4 ip) )) 
= tr ( E ((4 Il) ®4 n) )---(4 4p) ®4 4p) ))) 

iir" i*p 

3 ij 

= Y,^{ Ai r r)A( r r)n ) ■ 

In the case of periodic boundary conditions the SVD-based normalization can 
only be performed for all up to one matrix pair Ar and so the denominator 
matrix N r is non-trivial. However, the normalization is utilized to make the 
problem numerically stable (see |45|). 
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The matrix H r for the eigenvalue problem is given by 
E tr(4 4l) • • -A^) ■ tr(A« } • • -A?) • (e^Qf'e,) • • • (e^Q^e ip ) 

i±,...,i p 

i' lt -,i' p 
k 

E • • ■ 4^) ® (4° • ■ • a?)) • (g« >4l • • • g« ,< 

E tr [ E (^•* ) ®^ ) ))---(^-(^ ) ®4°); 

fc ii ,...,i p 

H,-,i' p 

(*0 C rffa) /l^ 



x H Hx = 



ll,...,t 

k 



EHnE(^-(4^^)) 



The effective Hamiltonian can be computed by contracting all indices except the 
indices representing the unknowns in matrix pair Ar , i r , i' r , m r , m r+ \ 1 k r , fc r +i 
leading to a 2D 2 x 2D 2 matrix H r . So the costs for solving the eigenvalue 
problem are in this case 0(D A ). In the periodic case also a SVD for N r has 
to be computed to solve the generalized eigenvalue problem numerically stable, 
which leads to costs of order D 6 . 



3.1.2. Matrix Product Density Operators 

The Matrix Product Operator approach extends the idea behind MPS from 
vectors to operators. Matrix Product Density Operators MPDO [46[ have the 
form 

Etr(4^ 4*'*>) W (22) 

with unit vectors e^. Matrix Product Operators MPO Q read 

E tr(4 ?l) A^^P h 9 — (8)P ip (23) 

with 2x2 matrices P, e.g. the Pauli matrices (JTJ) . Similarly, the tt format 
has also been extended to the matrix case (j47|)- These MPO concepts may be 
used for a representation of the Hamiltonian, such that the application of the 
Hamiltonian on an MPS vector leads to a sum of MPS vectors with less addends. 



1G 



3.2. 2D Systems: Approximation by Projected Entangled Pair States 

It is natural to extend the coupling topology of interest from linear chains 
to 2D arrays representing an entire lattice with open or periodic boundary con- 
ditions. To this end, the matrices in MPS are replaced by higher-order tensors 
thus giving rise to Projected Entangled Pair States (peps) [25|. Again, the un- 
derlying guideline is an area law (Zq . [HI . The extension to higher dimensions, 
however, comes at considerably higher computational cost: calculating expec- 
tation values becomes NP-hard (actually the complexity class is #P-complctc) 
[49I . This is one reason, why computing ground states in two-dimensional arrays 
remains a major challenge to numerics [5|. 

For 2D spin systems, the interaction between particles is also of 2D form, 
e.g., as described by Figure |U This leads to 2D generalization of MPS using a 




Figure 4: A 2D system in the open boundary case with physical indices ij and ancilla indices 
k and k. 

tensor network with small boxes related to tensor ( a^ lr ' B ^ ~ ~ ] with 5 legs. 
Thus, an inner product of two such vectors would look like 

t J k' r ,k r _^-^ ;k' s ,k s _^_ 1 k r ,k r +i ;k s ,k s +i 

In a PEPS vector, matrices are replaced by higher order tensors with one physical 
index related to the given vector x and the other ancilla indices related to nearest 
neighbors according to the lattice structure of the tensor network ([I(j). A peps 
vector can be formally written as 

x = Y.C R {{A I -})e lR , 

In 

where / stands for all a;-indices in Region R, Cr represents the contraction of 
indices following the nearest neighbor structure. 

To compute the inner product of two peps vectors, we have to find an efficient 
ordering of the summations. The related tensor network for the open boundary 



case is displayed in Figure 5(a) 
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In a first step all pairwise contractions relative to the vector indices i r>s are 
computed. Furthermore, in the produced boxes the indices are newly numbered 
by combining k rs with m r s to larger index k' r s . This generates Figure [5(b) | 



3- 



(a) Contracting physical indices. 



(b) Grouping index pairs. 



Figure 5: Contraction of two PEPS vectors: After the contraction of the physical indices related 
to the first column (a) the index pairs m rjS and k r . s are grouped to larger indices k' r s . 



Now the first and second column are contracted starting e.g. from the bot- 

Unfortunately, the boxes 



torn, resulting in the network displayed in Figure 6 (a 



in the newly generated left column have more indices than in the starting col- 
umn. So wc cannot proceed in this way directly. In order to keep the number 
of indices constant, an approximation step is introduced. The main idea is to 
reduce the number of indices by considering the first left column as MPS — with 
indices related to connections with the right neighbors as original physical in- 
dices (longer than 2) — and approximating the boxes by little tensors with only 
one leg instead of two to the vertical neighbors. 

Such a reduction can be derived by the following procedure. We want to 
reduce the rank of size D 2 in index pair {&2,ij ^2,2} to a new index k' 2 ± of 
size D. We can rewrite the whole summation in three parts, where c contains 
the contracted summation over all indices that are not involved in the actual 
reduction process. This leads to the sum 



E 



^{^2,3 '^2,1 }'{^3,1 '^3,2 } ^{^3,1 '^3, 



}>{*4 



2} ^1^3,3 '^4,1 '^42 }'f ^2,3 '^2.1 1 



The entries 



I (^2 



) 

'^2.2 ) ' (^3,1 '^'3,2 ) 



Collecting these matrices in a large rectangular matrix, we can apply the SVD 
on this matrix. Now we can truncate the diagonal part of the SVD to reduce the 
matrices A from size D 2 x D 2 to size D 2 x D. If we repeat this process along the 
first column all indices are reduced to length D. Note that this approximation 
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step is not affected by any accuracy demand but it is performed to keep the 
number of indices constant, which allows to proceed in an iterative way. 

Reducing the rank of the indices in the first column leads to the network 
illustrated by Figure |6(b)| with the same structure as Figure |5(b)| so we can 
contract column by column until we are left with only one column that can be 
contracted following Figure [7] 



[ 



(a) Contracting index pairs. 



(b) Result of the rank reduction. 



Figure 6: Contraction scheme for PEPS vectors: after the contracting illustrated in Fig. |5(b)| 
The newly generated first column in (a) has sets of indices leading to a more complex con- 
traction process as illustrated by the double lines. After the rank reduction (b), we are in the 
same situation as in Fig. |5(b)| and can thus proceed in the same manner. 




Figure 7: After application of the column-by-column contraction scheme, we end up with the 
last column, which can be contracted in the illustrated way. 

The overall costs depend on D 10 , resp. D 20 in the periodic case. The so- 
called virtual dimension D is usually chosen a priori by the physicists (e.g. 
D < 5), such that the computations can still be afforded ([45|). 
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|J3,1 

a 3,ly' 3 ,i j 2 ,l J2,2 



J2.1 



a 2.1;j 2 ,i,ji.i,ji,2 



a 2,2;j 2 ,2,ji,3ji,4 



Jl.l 



71.2 



°-l-2:j 



Jl,3 



J 1.4 



ai,4;j li4 ,i 7 ,<8 



Figure 8: Scheme of the Tree Tensor States (tts). 



3.3. Tree Tensor States and Multi-scale Entanglement Renormalization Ansatz 
In MPS the new spin particles are added one by one leading to a new matrix 
in the tensor ansatz. A modified procedure leads to Tree Tensor States (tts). 
In this ansatz the particles are grouped pairwise leading to a smaller subspacc 
of the original Hilbert space. This procedure is repeated with the blockpairs 
until only one block is left. This is displayed in Figure [5] and formula 

y ] (°i,iy'i,i.«i.«2 ' ' ' °Myi,4,»T,»8 ) (°2,iy'2,i Ji.i ,31,2 °2,2;j2,2,ji.3 j'1,4 ) (°3,iy3,i ) 
,32,2 

where each block is isometric: 

/~] a>k',i,j ■ a ks,j = 5k', k ■ 

id 

Hence, in this scheme the network consists in a binary tree built by small isomet- 
ric tensors with three indices. It can be shown that all MPS can be represented 
by tts (see [51|). 

A further generalization of tts and MPS is given by the Multi-Scale En- 
tanglement Renormalization Ansatz (mera, see [28|). Besides the top tensor t 
with two indices, the network is built by two types of smaller tensors: three 
leg isometric tensors (isometries) and four leg unitary tensors (disentanglers) . 
This formalism is displayed by Figure [9] In connection with eigenvalue approx- 
imation for mera, ALS cannot be used in order to optimize over the isometries 
that represent the degrees of freedom; instead other optimization methods have 
to be applied. 



4. Representations of Tensors 

As any state in a physical system (see Sec. [2]) with p particles may be seen 
as a p-th order tensor, we want to present some basics about tensors within this 
section. In the next section, we will give some modifications and generalizations. 
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18 '9 '10 '11 '12 '13 '14 «15 '16 



Figure 9: The Multi-scale Entanglement Rcnormalization Ansatz (mera) with periodic bound- 
ary conditions, consisting of isometrics ui and unitaries u. 



4-1. Some Basics about Tensors 

In its most general form, a pth-order tensor A = (a^,...^ ) G M. niX "' xn p 
is a multiway array with p indices. A first-order tensor is thus a vector, a 
second-order tensor corresponds to a matrix. If x and y are vectors (i.e. first- 
order tensors) it is well known that the outer product xoy := xy T is a rank-one 
matrix. As generalization, if oi, . . . ,a p are vectors, the tensor A :— a±o- ■ -oa pi 
which is defined as a^,...^ = '/; : , 02.,-. ■ ■ ■ a p; i p is a rank-one tensor. Hence, the 
application of outer products constructs tensors of higher order. The Kronecker 
product of matrices just corresponds to this definition of the outer product, but 
it reshapes the single tensor entries into a matrix of larger size. 

To begin, we write a tensor A as a sum of rank-one tensors: 



R 

A = Y,a[ j) oa { 2 j) o 



(24) 



If R is minimal in the representation (|24[) of A, we define the tensor rank of A 
as equal to R. 

For illustrating tensors, matrices are often in use. One possibility to bring 
back a general tensor to a matrix is given by the mode-n-unfolding (see (EH): 
When applying this technique, the tensor A = (a^....^^....^) is represented by 
the matrix 



A 



(n) 



-l» l n + l! 



,} ■ 
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The n-mode tensor matrix product is given as follows: Let A = (fflii,...,»„, ...,»„) 
be a p-th order tensor and U = a matrix, then the mode-n-product x„ 

is defined as 



A X n U — I S ' dii t ...,i n Uj t i 



Beside the total rank of a tensor there also exist a local rank concept: the 
mode-n-rank of a tensor is the rank of the collection of all vectors belonging to 
index n. If a given tensor A has the n-mode ranks r n , we define the multilinear 
rank of A as (r\ , . . . , r p ). 

4-2. Decomposition of Tensors 

For approximating tensors (xi lt ... t i ) there are two basic methods (see [52|): 
the Tucker decomposition and the canonical decomposition. The Tucker decom- 
position (HU) 

D 

•Eii,...,i p ^ ^ ymi,...,m p ^l;ii,mi ^2;i2 ; wi2 ^p\i p .m p (^5) 

mi , . ...m p 

represents the given tensor by a tensor y mi ,...,m p with less dimension in each 
direction; D is called the Tucker rank, y mi ,...,m p is called core tensor. This 



concept is illustrated in Figure 10(a) In the case of a binary tensor x £ R 2x '" x2 , 
this is not meaningful, because then the Tucker rank is already 2. 

The canonical decomposition (candecomp), which is also known as Parallel 
Factorization(PARAFAC), has the form 

^,..,* P =x>ii4i <4i- (26) 

s=l 

Hence, the PARAFAC decomposes the given tensor into a sum of rank one tensors 
(see [3 which is illustrated by Figure |lO(b)| If we think of a; as a vector 



this is equivalent to 

D 



with tensor products of smaller vectors. One often finds the normalized form 

D 



= X sa[ s) ®---®a^ (27) 



x 

s=l 



with vectors of norm 1. 

If D is minimal in the representation (|26[) of x, it is called the canonical 
rank. 
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Tensor Train Schemes 

Application of the concepts ([25]) or (|26| would be a generalization of SVD 
that allows a substantial reduction in number of parameters for deriving a good 
approximation for a given tensor. Unfortunately, these decompositions have 
disadvantages like still exponential growth, lack of robust algorithms for rank 
reduction. Oseledets and Tyrtyshnikov ([56[) proposed the following Tensor 
Train scheme (tt) as an attempt to overcome these difficulties. In a first step 
a dyadic decomposition is introduced by cutting the index set and the tensor in 
two parts introducing an additional summation with a newly introduced ancilla 
index to: 

x ii,...,ik 4 ,ik+i,--- t ip = ^~!°l;n,— ,«*."» ' 2;i*+i ip,m ■ (28) 



This is the first step of the Tree Tucker [55| decomposition. Now we may apply 
this process recursively to a\ and a?.. If we always use the index partitioning 
ij', ij+ii ■ ■ ■ i *pi we arrive at the tt format 



E 



m 2 ,...,iTip 



(29) 



Again we can distinguish between given, physical indices ij and ancilla indices. 



Figure 10(c) illustrates the tt decomposition concept. 
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(a) The Tucker decomposition scheme. 
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(b) The canonical decomposition scheme. 





m 2 


a 2;z 2 ,m2,m3 


m 3 


^3;'(3,m3,m4 


m 4 








1 

h 




1 

«2 




1 





(c) The Tensor Train decomposition scheme 
Figure 10: Tensor decomposition schemes 
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The tt scheme (f2l)|) is exactly the MPS form ^ for open boundary condi- 
tions: 

Xi u ... )ip =G^G^ (30) 

with matrices Gj of size Dj x -Dj+i, (D\ = D. p+ i = 1), where the matrix sizes 
differ and are called the compression rank. 

In (5^ [. Khoromskij generalizes the tt concept to Tensor Chains via the 
definition 

•Eil,...,i p ^ ^ ,mi ,m2 ^2;i2 ,^2 ,«i3 

mi ,...,m p 

which corresponds to MPS with periodic boundary conditions, i.e. 

x h ,..., ip =tr(G^G^ GM) . (32) 

Starting with Formula (|28[) this step can also be applied recursively to the 
two newly introduced tensors (04) and (a 2 ) with any cutting point. Cutting a 
tensor in two parts should introduce as many ancilla indices as described by the 
tensor network connecting the two parts. So in ID there is only one neighboring 
connection that is broken up and introduces 1 additional index. 

In the 2D network described in Figure |4] wc could apply this method by 
cutting in a first step the down left knot in, thus introducing two ancilla indices 
labeled with (21,1,21,2), resp. 1,12,1), representing the broken connections to 
the neighbors of ii 1. Proceeding in this way knot by knot leads to the peps 
form. But generalizing the approach we can also consider general cuts of region 
R of the tensor network in two regions R\ and i?2, replacing the given tensor 
by a sum over a tensor product of two smaller tensors with indices related to 
i?i, resp. i?2, introducing as many ancilla indices as broken connections in the 
cut. 



*p;a p ,m p ,mi 5 



(31) 



5. Modifications of parafac and MPS Representations 

In this section we want to introduce some own ideas concerning the general- 
ization and modification of both the PARAFAC decomposition scheme (see 14. 2|) 
and the MPS (tt) format. These newly presented formats will be used as vector 
ansatz for the minimization of the Rayleigh quotient ^ . As we have pointed 
out in Section 12. 3[ we are looking for representation schemes that allow both 
proper approximation properties and efficient computations such as fast con- 
tractions for inner product calculations and cheap matrix-vector evaluations. 
In view of these requirements we present both theoretical results concerning 
the computational complexity and numerical results showing the benefit of our 
concepts. 

At this point we want to emphasize that our problem does not consist in 
the decomposition of a known tensor but to find an appropriate ansatz for 
vectors and to formulate algorithms for the ground state computation working 
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on such vector representations. Hence, in this context we focus on algorithmic 
considerations and investigate, which modifications and generalizations of the 
presented concepts are still affordable to solve physical problems, where we a 
priori work with approximations with ranks which are given by the physicists. 
It turns out that the usage of low-rank approaches suffices to give proper results, 
compare, e.g., [BH ]. 

5.1. PARAFAC Formats 

Any state x € C 2P of a physical system can be tensorized artificially in 
several ways. In the easiest way we may rewrite a; as a pth ordered binary 
tensor of the form 



But we may also define blockings of larger size, which will follow the interac- 
tions of the physical system (e.g. nearest-neighbor interaction). These blocking 
concepts introduce formats with a larger number of degree of freedoms promis- 
ing better approximation properties but still allow efficient computations of 
matrix-vector and inner products. For reproducing the physical structure such 
as internal and external interactions we will also allow formats with different 
overlapping blockings in the addends, compare Subsection 15.31 

Such blockings of indices can be seen as tensorizations of lower order q < p: 



where the fcth mode combines tk := — s^—i binary indices and has therefore 
index jj, = (i Sk _ 1 +i, . . . ,i Sk ) of size 2 th (for reasons of consistency we define 
so = and s q = p) . 

One way to find a convenient representation is to consider appropriate de- 
compositions of the tensorization (f34|) . In this context we don't consider the 
Tucker format (|25p as it is only meaningful to decompose very large mode sizes 
(meaning large sets of blocked indices). 

Hence, from now on we consider the PARAFAC concept and choose the ansatz 
vector to be a sum of rank- 1 tensors (see Figure ITTj) . In the simplest case (f3"3"|) , 
the PARAFAC decomposition takes the form 



a sum of tensor products of length 2- vectors. In view of (fTTj) , the decomposi- 
tion ([33)) can be seen as a special MPS form. Indeed, every PARAFAC represen- 
tation (|3"5"|) with D addends corresponds to an MPS term with D x D diagonal 
matrices. This fact becomes clear from the construction 



«E i-^ii ... .,ip^)ia— 0.1 • 



(33) 



x — { x (i 1 ,...,i 31 ),...,(i 3q _ 1+1 ,...,i Sq )) — { x ji,....j q ) ■ 



(34) 




(35) 





A (o) 

r 
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More generally, the PARAFAC scheme for the tensorization leads to the 
ansatz 

D 

x = Y J x i ) ®---® x f ■ ( 36 ) 

with vectors x^p of moderate size 2 tj . Figure [11] illustrates such a decomposi- 
tion. 



X 




Figure 11: PARAFAC ansatz for a chosen tensorization of the vector x to be represented. 



Computational Costs 

Inner product calculations of two representations of the form (|36p reduce to 
the inner product of the block vectors: 

y*x = £ (y^xf'U^xP) • ■ • (yW H xf )) . (37) 

Therefore, the total costs for each of the D 2 inner products are 

2(2*! + 2* 2 + • • • + 2 t ") + q < q(2l + 1) < p(2l + 1) 

where t := maxjii, t2, ■ ■ ■ ,t q } is the maximum number of grouped indices and 
thus I = 2* denotes the largest block size in a;. If all q index sets have the same 
size t (i.e. t = p/q), the costs can be bounded by q(2 ■ T>li + 1). 

To compute the matrix-vector product efficiently, we group the Hamilto- 
nian H ^ in the same way, i.e. 

M 

h = Y, MQ { i k) ® • ■ • ® ® • • • ® ® • • • ® Qp fc) ) 
fe=i 

M 
fc=l 



2G 



For the matrix- vector product we thus obtain 

D 

Hx 




M D 



fc=i e=i 

M D 



= EEwi M) «-®»i M> - 

k=l 1=1 

However, typically we do not need the matrix-vector products psp explicitly, 
we only require them for inner products of the form y H (i?x) (compare the 
nominator of the Rayleigh quotient) 

m d / „ x 

fe=i i,i'=i \ ' 
The products H^x^ can be computed implicitly without constructing the 

(k) 

matrices H- explicitly. Each of these small matrix-vector products can be 
computed linearly in the size of (i.e. 2 4< ). Thus, the total costs for each 
addend in the inner product (|39|) are 3(2 41 + 2* 2 + • • ■ + 2 tq ) + q and can again 
be bounded by q(3 • 2 p / q + 1) in the case of equal block sizes. Hence, the costs 
for (J3SJ) are 2MD 2 q{3 ■ 2 p li + 1). 

Using PARAFAC for the Ground State Problem 

Let us now apply the PARAFAC approach (|3"6")) as ansatz for the Rayleigh 
quotient minimization. Then, Eq. [6] reads 

H 



£ ® •■ • ® xMj H (f: xf } ® ■ •■ ® xWj 



(40) 



This minimization task can be realized by an ALS-based procedure. 

In a first way we could think about the following proceeding: We start with 
a PARAFAC representation of rank 1 (D = 1), optimize it via ALS and then we 
successively add one summand and optimize it in an ALS-based way. In the first 
stage we would have to minimize 

. x H Hx . (xi ® x 2 ® • • ■ ® x„) H H (xi ® x 2 ® • • ■ ® x„) 

mm — 77 — = mm 



a: X"X g, ^ g . . . ( Xl ^2 . . . (g, aj^) 

This optimization problem can be solved via ALS: considering all of the a)j up 
to Xi as fixed, we obtain 

. (Xx ® ■ ■ ■ ® Xi (g) ■ ■ ■ ® Xa)^ H (Xx ® ■ ■ ■ ® Xi <E) ■ ■ ■ ® x„) 

mm — 7h~, — ■ ( 41 ) 

Xi [X\ ® ■ • ■ ® Xi <g> • • • ® Xq) (Xl ® • • • ® Xi ® • ■ • ® Xq) 
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Following and (|57)) . we may contract all indices up to i and obtain 

E a fc (*i H #f W) • ■ • (xi H fff ^i) • • • (x q H H^x v ) 
xi (a;i H cci) • • • (sc^a;*) • • • (x^) 

M , . / M 

E akpkxVH^Xi Xi B E ^iff' 

fc=i • Vfc=i 



a standard eigenvalue problem for a matrix of size 2 fi x 2 4i that can be solved via 
classical iterative methods which only require the computation of matrix vector 
products. These products can be executed implicitly without constructing the 
matrices explicitly. 

Now we suppose that we have already optimized D — 1 addends in the 
representation (|36[) and we want to find the next optimal addend. This means 
that the ansatz vector is now of the form 



x = x 1 ®---®x q +Y~' y[ £) <g) ■ • • <g> y<p 



= (JS> <*>j 

3=1 



with already optimized y y - -terms and vectors xi that have to be optimized via 
ALS. Contracting over all terms up to Xi (compare Eq. we obtain 

x A Hx = Xi H Hi 

where Ui and /3 comprise the contractions with the yj terms. For the denomi- 
nator we analogously obtain 

X K X = Xi^jIXi + Vi H Xi + Xi^Vi + p . 

Altogether we have to solve the generalized eigenvalue problem 



Hi Ui \ ( Xi 

p ) I i 



X* 1 



■yl Vi \ ( Xi 
Vi* p { 1 



(42) 



It turns out that the denominator matrix can be factorized via a Cholesky 
factorization. Therefore we have to solve a standard eigenvalue problem of 
moderate size. 

However, the proposed procedure may cause problems which result from the 
fact that, for general tensors, the best rank-D approximation does not have to 
comprise the best rank-(Z? — 1) approximation, see, e.g., 52j. Our numerical 
results (Figure fT2|) will approve this fact. 
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For this reason we now consider a technique where we use ALS to optimize 
all the blocks related to the same mode i simultaneously. The minimization 
task takes the form 



D 



E x i 



X 



(A) 



x 



(E> 



x 



(A) 



x 



E4« 



.(<) 



D 

E^ J 



.(*) 



Contracting all terms in the nominator (compare Eq. l39| up to the x\ vectors 
results in 



M D 



k=i i,i'=i 

D 



E« fc E (*P *l h) *? 



(A') rrW W 

x\ H\ >x\ 



E «!° H W 



For the denominator, we analogously obtain 

D 



E *f )H (/w)* 



9 q 



Altogether, the approach leads to the minimization problem 



I fr(f ,1) 



(43) 



a generalized eigenvalue problem of size D2 ti x Z)2*\ This approach requires 
the solution of larger eigcnproblems, but by increasing the set of variables to be 
optimized it turns out that we require less optimization steps and obtain better 
convergence results (see Figure [THJ) . 



Formulation of an Algorithm 

As possibility to choose the initial guesses, we propose to set {x^ , . . . , x\ D ^ } 
to D linearly independent eigenvectors of 

M 



Hi = E Qs^+i 



® ...<gll 



)(fc) 



k=l 
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Algorithm 1 Raylcigh Quotient Minimization in the parafac model 

1: Provide initial guess for all vectors x\ (i = 1, . . . , q and £ = 1, . . . , D) 
2: while Not yet converged do 
3: for i = 1, . . . , q do 

4: Compute all contractions in x K Hx and x^x up to index i 
5: Solve the generalized eigenvalue problem (j43|) 
6: end for 

7: Normalize each addend in the PARAFAC representation (|27|) 
8: end while 



As stopping criterion we may define a maximum number of iterations or we 
could specify a stopping criterion based on the degree of improvement. One 
possibility to accelerate the ALS-based optimization procedure is to apply the 
Enhanced Linesearch method [H^|, which is not considered here. 

Numerical Results 

After these technical derivations, we want to show first numerical results. We 
computed the ground state of a 10 spin and a 12 spin Ising-type Hamiltonian ([3]) 
using different blocking schemes [t\, . . . ,t q ]. The ALS-based Rayleigh quotient 
minimization was performed in two different ways: a one-by-one minimization 
following (|42[) and a simultaneous minimization described by Algorithm [T] It 
turns out that the results are getting better when using larger vector blocks. 
This consideration is based on the fact that using larger blocks leads to a larger 
number of degrees of freedom. Figure [12] depicts the effect of using different 
PARAFAC representations and shows the benefit resulting from the simultaneous 
optimization. 

5.2. Block MPS 

In the block MPS ansatz we group single particles together and apply MPS 
concepts to these grouped blocks. The MPS ansatz (j9|) is based on using 2 
matrices A^ 3 ^ at position j. A first generalization of MPS vectors can be derived 
by blocking indices in subgroups. So we can combine indices i r ) to index 

ji, and indices (i r +i> ■■■tip) to index Then the MPS vector to this blocking 
with open boundary conditions would give 

T — T- ■ — Tr , r 4^ . 

i — J >ii,...,i p — ■ i {ii,...,i r },{i r+ i,...,i p } — ^31,32 — -"-1 ^2 

or vector 

D 

x = ^ ai ; i, m2 ® a 2; m 2 ,i (44) 

7712 — 1 

with D pairs of short vectors of length 2 r , resp. 2 p ~ r . Compared to standard 
MPS we have different numbers for the degree of freedom (DOF), and different 
costs for inner products. MPS combines in each term p vectors of length 2 in a 



30 




- blocking [3, 4, 3] 
-blocking [4, 2, 4] 

- blocking [5, 5] 

- blocking [G, 4] 

2 3 4 

Number of addends 




-blocking [3, 4, 3] 
-blocking [4, 2, 4] 

- blocking [5, 5] 

- blocking [G, 4] 



3 4 
Number of addends 



(a) 10 spins, one-by-one optimization. (b) 10 spins, simultaneous optimization. 

10" i 




a 10- 



- blocking [3, 3, 3, 3] 

- blocking [3, G, 3] 
-blocking [4, 4, 4] 

- blocking [6, 6] 




- blockin; 

- blocking [3, 6, 
-blocking [4, 4, 

- blocking [0, G] 



2 3 4 

Number of addends 



3 4 
Number of addends 



(c) 12 spins, onc-by-one optimization. (d) 12 spins, simultaneous optimization. 

Figure 12: Approximation error for the computation of the ground state energy of a 10 spin 
and 12 spin Ising-type Hamiltonian. 



tensor product and combines D p vectors; thus, it has 2pD 2 DOF's and 2pD 3 
costs for inner products. The 2-term block coarse grain MPS form (|44|) combines 
in each term 2 vectors of total length 2 r and 2 p ~ r . Therefore, it has D(2 r + 2 p - r ) 
DOF's and the estimate for the inner product costs is D 2 (2 r + 2 p ~ r ). The 
computation of the full inner product reduces to D 2 inner products of short 
vectors at position 1, resp. 2: 

D D 

X^X = ^2 ai;l,m' 2 H ® a2;m4,l H ' ^ Ol;l,m 2 ® a 2 ;m 2 ,l 

m' 2 «t2 

= ( a l;l,m' 2 H a l;l,m 2 ) ' (a2;m^,l Ha 2;Tn 2 ,l) • 

m,2 ,m 2 

For using three blocks we derive 

T- — Xr , r t r Ah) Ah) Ah) 

■Li — J '{li,...,l ri },{l ri + i,...,l r2 },{*r 2 + l,...lp} — X 3l,32,33 ~ ^1 ^2 ^3 
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or vector 

D 

X = ai;l,m a ® a 2;m,2,m 3 &> fl3;m 3 ,l 

m2 ,m3 

with D2 fl + D 2 2 t2 + D2 t3 DOF's, h = n, t 2 = r 2 - n, t 3 = p - r 2 , combining 
D 2 long vectors, and costs of order D 4 for an inner product in view of 

X X = { a l;l,m' 3 0-l;l,m 2 ) ' ( a 2;m' 2 ,m' 3 a 2;m2,m 3 ) ' { a 3;m' 3 ,l a 3,rra 3 ) ■ 



By using the MPS contraction scheme (Figures 2(a)| — |3(d)| ), the costs for the 



inner product can be reduced to D 3 (2 tl + 2* 2 + W) 
In general, the Block MPS ansatz with k blocks gives 

T- — T- — T- ■ — . 4O2) . , . 40'*-l) . 4O'*) (AZ.\ 

%l — ,J„ — x ]i,... ,3k — A l A 2 A k-1 A k y^ ) 

or vector 

X = ^2 ai;l, ma ® a2;m 2 ,m 3 ® ■ • ■ a k -l;m k _ 1 ,m h ® afc;m fc ,l ■ (46) 
m2,-"" 5 m fc 

DOF's: L»(2 tl + 2**) + L> 2 (2* 2 + • • • + 2**- 1 ) for ti = n, t fc = r fe - r fc _i. 
Costs: D 3 (2* 1 + • • • + 2 tk ) using the MPS contraction scheme, 
Combination of D k ~ 1 full vectors. 

Altogether, the Block MPS ansatz causes greater effort for the computation 
of inner products, but the increasing number of degrees of freedom may allow 
better approximation properties. The Block MPS format corresponds to the 
Tensor Train format with mode sizes 2 li , i = 1, . . . , k. 

In summary, MPS can be considered as a fine grain method, using many 
vectors, while the coarse grain MPS uses only a sum of tensor products of a 
few large vectors respectively. The fine grain MPS uses at each position s a 
matrix pair As* 3 ', i s = 0, 1, while the fc-form (|45p uses 2* s matrices As , j s = 
0, 1, 2 ts — 1. In the same way the length of vector a S ;m s ,rri s+1 is 2* s . 

5.3. Mixed Blockings in the PARAFAC Ansatz 

In this section, we want to analyze mixed tensor representations similar to 
PARAFAC concepts presented in Section 15.11 but allowing different blockings in 
the addends (see Figure I13[) . This may also be interesting for mixing, e.g., 
Block-MPS vectors with different block structure (compare Subsection 15. 4|) . 

Reasoning for mixed blockings: In connection with PARAFAC or MPS, 
approximations of a given vector as a sum of tensor products of smaller vectors 
(such as (flTj) or ([50])) which are of different sizes might give different approx- 
imation qualities. E.g. using a partitioning with two vectors is related to a 
certain matricization of the given vector. The matricization leads to an SVD 
with related singular values of certain decay behavior and therefore a certain 
truncation error depending on the number of singular values that are used for 
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Figure 13: A PARAFAC-like representation. In contrast to the decomposition scheme illustrated 
by Figure [TT1 different blockings are possible. 

the approximation. Hence, it might be useful to combine different matricizations 
in PARAFAC allowing tensor products of different partitionings. 

In the following we want to discuss the question which blockings allow effi- 
cient contractions similar to the case of uniform blockings. Using such mixed 
schemes for finding ground states also leads to generalized eigenproblcms sim- 
ilar to (|42I43[) . but now the computation of inner products is generally more 
costly. In contrast, the application of the Hamiltonian ([2]) to a mixed vector 
representation can again be computed efficiently because the Hamiltonians are 
given by Kronecker products of 2 x 2 matrices. The question here is, what type 
of blocking allows fast and efficient contractions in inner products for vectors 
with different tensor blocking. 

In the following we will consider the problem what tensor aware permutations 
can be combined that allow fast and efficient contractions. Such permutations 
allow a generalization of PARAFAC, where in each approximation step a new 
tensor aware permutation can be applied before the new term is computed. 

So in the first, elementary case we consider 

z = {z iu ...,i p ) = (xi. tilt ... tiri )<»••• <g> (a;fc;* rk _ 1+1 ,...,* p ) 

+ (yi;h,...,i B1 ) ® ••• ® (2/m;», m _ 1+1 ,...,i„) 

where the vector z with 2 P components is decomposed in two vectors that 
are Kronecker products of smaller vectors. The length of the k vectors Xj = 
( a; j)*r-_i+i,.",*r-) is given by 2 rj_r ^ 1 and similarly the m vectors yj areoflength 
2 s i~ s '- 1 . The blocking of x and y may be different. Such forms can be seen as 
generalizations of PARAFAC, or the coarse grain MPS form. 

The interesting question is whether it is possible to introduce efficient con- 
traction schemes for such sums. So in the following we want to analyze the 
costs for inner products of two vectors that are built as tensor products, but 
with different block structure as displayed in Figures UM and [15l 
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The One- Dimensional Case 

If both vectors have the same block structure, the inner product reduces to 
the inner product of the block vectors, compare (|37|) . The costs can then be 
bounded by 3fc2P/ fc . 

The more interesting case appears if we allow different block sizes. First, let 
us consider the open boundary case where the first block in both vectors starts 
with i\ and the last block ends with i p , compare Figure [T4l The contraction 



Vi+i 



















'S m -l + l *s m 



Figure 14: Contraction scheme for the calculation of the inner product of two vectors, which 
are formed by tensor products of smaller vectors, in the open boundary case. 



starts with the left block computing the summation over i\, ... 



E 



2/1; 



(48) 



This leads to an update of y\ which is replaced by a new vector to indices 
Vi+ii •••) *«!• The costs are 2 ri for one summation, and the summation has to 
be done for all indices i ri +i, ■■-,i Sl , hence 2 Sl ~ ri times, leading to total costs of 
2 Sl for this step. Now we have to repeat this procedure for the pair X2-i ri+1 ,...,i r2 
and the updated • Let us assume that all blocks have size bounded 

by 2 r . In each step the costs are bounded by the size of the longest block, and 
the block size of the intermediate generated vectors is smaller than the size of 
the previous vectors. Hence in total the costs are bounded by the number of 
necessary contractions = k + m times the maximum block size: 2 r ■ (k + m). 

In the periodic case, the first and the last blocks are allowed to overlap as 
displayed in Figure [T5] In this case the first contraction looks like 



Figure 15: Contraction scheme in the periodic boundary case. Differing to open boundary 
case illustrated in Figure [Ml the first and last blocks also overlap. 



^ ' 2'X;ti,...,i» 1 ,...»r 1 ' ym;is m + i,...,ip,ii,...,is 1 ^ • 2 l;t, m +i ipAn + i *rj • (49) 

ii •••»»! 

The costs are given by the length of the contraction times the number of times it 
has to be executed, that is 2 Sl • 2 p ~ Sm • 2 n ~ Sl . In the worst case the summation 
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part could be only one index which leads to 2 2r 1 costs for maximum block size 
2 r , where r := max{ri,r 2 - r\, . . . ,r k - r fc _i,si,s 2 - s X) . . . ,s m - s m _i}. But 
in each step we can choose a contraction where the summation part is larger 
than the remaining indices. Then, the costs for one contraction are bounded by 
2 ri+( P - Sm ) < 2 3r/2^ Therefore, in total the costs are bounded by 2 3r / 2 (m + k). 

5.4- Mixed Blockings for Block MPS Vectors 

The contraction schemes (|48| presented in the previous subsection may also 
be applied to calculate inner products of Block MPS vectors (|45|) as introduced 
in Section 15.21 To this end we have to append the ancilla indices resulting from 
the matrix products to the physical indices. For contracting the inner product 
of the Block MPS vectors 

x = £ 4 

»1,...,»J> 

y= E *? 

ii,...,i p 

we can proceed in the same way as for standard MPS vectors (see Figures 13]) , 
but instead of contracting over only one physical index ij per step, we have to 
sum over all overlapping indices as illustrated in Figure 1161 
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■ >*«!, »«! + !, ■ 


. , i ri , hi, A/2 
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.,i T2 ,k 2 ,k 3 





...,i ri )^(ir 1 + i,...,V 2 ) _ _ _ ^(ir- k _ 1 + i,---,ip) 



(50) 



■ • • £?„ 




Figure 16: Contraction scheme for the inner product of the block MPS vectors 1501 . 



5. 5. The parafac Model for the Two- Dimensional Case 

For a 2D tensor network, we introduce another block structure that also 
allows fast contractions. As example we consider a 2D physical system with 
nearest-neighbor interaction in both directions as displayed in Figure 1171 

In a first step we group r indices together to larger subblocks taking into 
account the neighborhood relations, see Figure [l8(a)[ 

We denote the newly introduced blocks by indices j 8 . In the above example 
each j block contains 4 physical i indices (Figure ["l8(b)[ ). 

Based on the introduced coarse numbering, we define a final collection of 
blocks by building pairs of subblocks. We allow four different patterns of su- 
perblocks (see Figure H~9| , which enable efficient contractions in the calculation 
of mixed inner products. 
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Figure 17: A 2D physical system with nearest neighbor interaction in both directions and 
periodic boundary conditions. 



These four patterns are related to four different vector types 



x l 


= a jl,32 ^ 


' a js,ji 


® a j5 , J6 




® a jl3,jl4 & 


a jl5,jl6 
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= a j4,jl $ 


^ a jz,3s 


® a is ,3s 




® a jl6,jl3 & 


a jl4,J15 


x :i 


= a 3l,3s @ 


^ a 32,je 


® a 33,3r 


® • • 




a jl2 ,J18 


x l 


= a jl3,jl 


^ a jl4,j2 ® a jl5 


,J3 ® 


■•■(g) Oj- T 


® a j8,jl2 



(51) 



As ansatz vector we may consider combinations of the four vector patterns: 
2 = ^ Qt,as; + ^ fta^ + ^ + ^ <J K a£ . 

Contractions between vectors of type 1 — 4 can be done in C(2 3r ) operations. 
Each contraction considers a block of vector x s combined with a block of vector 
Xt that has an index jk in common. Therefore, the costs are the subblocksize 
to the power of three. The result of such a contraction is a new block consisting 
of the pair of indices that were left over. Hence the total costs are the total 
numbers of blocks times the subblocksize to the power of three. 

Remark: Similar to the generalization of MPS to Block MPS with mixed 
blockings as introduced in Section [5~4l the mixed blocking formats for 2D sys- 
tems could also be applied to modify the peps ansatz (see Section l3~2"T) . For effi- 
cient computations of inner products we may generalize the contraction scheme 
(Figures [5] and |6|) to sets of overlapping physical indices. 



6. Conclusions 



This work combines similar concepts proposed in quantum information and 
linear algebra communities and presents them from a (multi-) linear algebraic 



3G 




(a) Grouping of indices. 



(b) Introduction of new indices. 



Figure 18: The 2D physical system as introduced in Figure [T71 but now adjacent indices are 
grouped together and denoted by new indices j. The original neighborhood relations from 
Figure fT7l arc still taken into account. 
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Figure 19: Four different patterns of subblocks related to the four vector types defined in 15U 



point of view in a unifying mathematical and computational framework. In this 
context we are mainly interested in algorithmic considerations. 

For our main task, the computation of ground states of physical quantum sys- 
tems, we have to minimize the Rayleigh quotient over an exponentially growing 
vector space. To overcome the curse of dimensionality we have to find appro- 
priate formats which allow cheap computations of matrix products and inner 
products, but still guarantee proper approximation properties. These formats 
are structured in such a way that they enable relaxation methods such as Al- 
ternating Least Squares. 

Physicists have developed concepts like Matrix Product States, which cor- 
respond to Tensor Trains, for linear one-dimensional problems or Projected 
Entangled Pair States for two-dimensional problems. We present these con- 
cepts from a mathematical point of view and show how computations such as 
contractions of inner products can be performed efficiently. 

As an ansatz format to minimize the Rayleigh quotient we also consider 
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the PARAFAC format, which allows modifications in several directions, such as 
different blocking structures. These generalizations, which can also be applied 
to modify, e.g., the MPS or PEPS ansatz, are constructed in such a way that 

f. the physical (inter) actions (e.g. nearest neighbor interactions) can be 
reproduced properly, 

2. representations are only polynomial in the system size p, 

3. computations such as inner product contractions can still be performed 
efficiently. 

Our mixed blocking ansatz in PARAFAC and MPS is an efficient and flexible 
way to represent states of physical systems and can easily be extended to higher 
dimensions or systems with more complex interactions. It is, however, an open 
question, how the difference in the blockings can be constructed automatically 
and how efficient they can be for general problems. 
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